References

Background

Here we have 82 data points in the galaxy from a 1-dimensional unknown distribution. The aim is to fit a normal finite mixture model with a pre-specified number of latent clusters.

Load packages

library(tidyverse)
library(rstan)
## devtools::install_github('jburos/biostan', build_vignettes = TRUE, dependencies = TRUE)
## library(biostan)
library(DPpackage)
library(tidybayes)
library(bayesplot)
set.seed(732565397)

Prepare data

data(galaxy, package = "DPpackage")
galaxy <- galaxy %>%
    as_data_frame() %>%
    mutate(log_speed = log(speed),
           k_speed = speed / 1000)
## Warning: `as_data_frame()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.
galaxy
## # A tibble: 82 x 3
##    speed log_speed k_speed
##    <dbl>     <dbl>   <dbl>
##  1  9172      9.12    9.17
##  2  9350      9.14    9.35
##  3  9483      9.16    9.48
##  4  9558      9.17    9.56
##  5  9775      9.19    9.78
##  6 10227      9.23   10.2 
##  7 10406      9.25   10.4 
##  8 16084      9.69   16.1 
##  9 16170      9.69   16.2 
## 10 18419      9.82   18.4 
## # … with 72 more rows
ggplot(data = galaxy, mapping = aes(x = k_speed)) +
    geom_point(y = 0.5) +
    scale_y_continuous(limits = c(0,1), breaks = NULL) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

The speed data seem to show some distinct clusters. We will use the following grid for later visualization.

grid_max <- 40
grid_min <- 0
grid_length <- 120

Define helper functions

Helper functions for here and later use.

print_relevant_pars <- function(fit, pars = c("mu","sigma","Pi","lp__")) {
    print(fit, pars = pars)
}

traceplot_all <- function(fit, pars = c("mu","sigma","Pi","lp__")) {
    for (par in pars) {
        print(traceplot(fit, inc_warmup = TRUE, pars = par))
    }
}

pairs_plot_all <- function(fit, pars = c("mu","sigma","Pi")) {
    for (par in pars) {
        pairs(fit, pars = par)
    }
}

plot_draws <- function(stan_fit, chains = NULL) {
    ## Note direct access to global variables
    draw_data  <- tidybayes::tidy_draws(stan_fit) %>%
        select(.chain, .iteration, .draw, starts_with("log_f")) %>%
        gather(key = key, value = value, starts_with("log_f")) %>%
        mutate(key = gsub("log_f|\\[|\\]", "", key) %>% as.integer(),
           x = factor(key, labels = seq(from = grid_min, to = grid_max, length.out = grid_length)) %>%
               as.character() %>%
               as.numeric(),
           value = exp(value))

    if(!is.null(chains)) {
        draw_data <- draw_data %>%
            filter(.chain %in% chains)
    }

    summary_density <- draw_data %>%
        group_by(.chain, x) %>%
        summarize(value = mean(value))

    ggplot(data = draw_data, mapping = aes(x = x, y = value,
           group = interaction(.chain, .iteration, .draw))) +
    ## geom_line(size = 0.1, alpha = 1/20) +
    geom_line(data = summary_density, mapping = aes(group = .chain), size = 0.5, color = "gray") +
    geom_point(data = galaxy, mapping = aes(x = k_speed, group = NULL), y = 0) +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())
}

Finite Mixture of Normals

Finite Mixture of Normals (Unordered mu)

3 latent cluster model

H <- 3
normal_fixed_mixture_stan_tau3 <-
    rstan::stan(model_code = readr::read_file("./bayesianideas_density_normal_fixed_mixture_tau_unordered.stan"),
                data = list(
                    ## tau ~ Gamma(alpha, beta)
                    alpha = rep(10^(-3), H), beta = rep(10^(-3), H),
                    ## mu ~ N(m, s_squared)
                    m = rep(20,H), s_squared = rep(1000,H),
                    ## Pi ~ Dirichlet(dirichlet_alpha / H)
                    dirichlet_alpha = rep(1,H) * 1 * 1,
                    ## Number of clusters
                    H = H,
                    n = nrow(galaxy),
                    y = galaxy$k_speed,
                    grid_max = grid_max,
                    grid_min = grid_min,
                    grid_length = grid_length),
                iter = 2000,
                chains = 12*3,
                seed = 34618604)
## Warning: There were 13 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Check lp__ traces to detect better fits
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::gather_draws(lp__) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Order fits by mean lp__ to find bad ones
normal_fixed_mixture_stan_tau3_good <-
    normal_fixed_mixture_stan_tau3 %>%
    tidybayes::gather_draws(lp__) %>%
    group_by(.chain) %>%
    summarize(mean = mean(.value)) %>%
    arrange(mean) %>%
    ## Define
    mutate(good_chain = between(mean, max(mean) - 1, max(mean)))
normal_fixed_mixture_stan_tau3_good %>%
    print(n = Inf)
## # A tibble: 36 x 3
##    .chain  mean good_chain
##     <int> <dbl> <lgl>     
##  1     15 -226. FALSE     
##  2     21 -213. FALSE     
##  3     17 -211. TRUE      
##  4     24 -211. TRUE      
##  5      2 -211. TRUE      
##  6     30 -211. TRUE      
##  7     36 -211. TRUE      
##  8      1 -211. TRUE      
##  9     16 -211. TRUE      
## 10     19 -211. TRUE      
## 11      7 -211. TRUE      
## 12     26 -211. TRUE      
## 13      3 -211. TRUE      
## 14     20 -211. TRUE      
## 15     10 -211. TRUE      
## 16     33 -211. TRUE      
## 17      6 -211. TRUE      
## 18      9 -211. TRUE      
## 19     25 -211. TRUE      
## 20     27 -211. TRUE      
## 21     22 -211. TRUE      
## 22     34 -211. TRUE      
## 23     12 -211. TRUE      
## 24      4 -211. TRUE      
## 25     14 -211. TRUE      
## 26     29 -211. TRUE      
## 27     35 -211. TRUE      
## 28     31 -211. TRUE      
## 29      8 -211. TRUE      
## 30      5 -211. TRUE      
## 31     13 -211. TRUE      
## 32     23 -211. TRUE      
## 33     32 -211. TRUE      
## 34     11 -210. TRUE      
## 35     28 -210. TRUE      
## 36     18 -210. TRUE
## For vague prior fit
bad_chains <- normal_fixed_mixture_stan_tau3_good %>%
    filter(!good_chain) %>%
    magrittr::extract2(".chain")

## Check lp__ traces for good chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::gather_draws(lp__) %>%
    filter(!(.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Check lp__ traces for bad chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::gather_draws(lp__) %>%
    filter(.chain %in% bad_chains) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Rhat for good chains (slightly different definition?)
lapply(normal_fixed_mixture_stan_tau3@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,-1 * bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 1.002016
## Rhat among bad chains
lapply(normal_fixed_mixture_stan_tau3@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 4.143634
## Estimated densities among good chains
plot_draws(normal_fixed_mixture_stan_tau3,
           chains = setdiff(seq_len(length(normal_fixed_mixture_stan_tau3@sim$samples)),
                            bad_chains))

plot_draws(normal_fixed_mixture_stan_tau3,
           chains = setdiff(seq_len(length(normal_fixed_mixture_stan_tau3@sim$samples)),
                            bad_chains)[1])

## Scatter plot of parameter values
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::tidy_draws() %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = `mu[1]`, y = `mu[2]`)) +
    geom_point() +
    facet_grid(~ bad_chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among good chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among bad chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among good chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among bad chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among good chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among bad chains
normal_fixed_mixture_stan_tau3 %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

pairs_plot_all(normal_fixed_mixture_stan_tau3)

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

pairs(normal_fixed_mixture_stan_tau3, pars = c("mu","sigma","Pi"))
## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

plot_draws(normal_fixed_mixture_stan_tau3)

## shinystan::launch_shinystan(normal_fixed_mixture_stan_tau3)

6 latent cluster model

H <- 6
normal_fixed_mixture_stan_tau6 <-
    rstan::stan(model_code = readr::read_file("./bayesianideas_density_normal_fixed_mixture_tau_unordered.stan"),
                data = list(
                    ## tau ~ Gamma(alpha, beta)
                    alpha = rep(10^(-3), H), beta = rep(10^(-3), H),
                    ## mu ~ N(m, s_squared)
                    m = rep(20,H), s_squared = rep(1000,H),
                    ## Pi ~ Dirichlet(dirichlet_alpha / H)
                    dirichlet_alpha = rep(1,H),
                    ## Number of clusters
                    H = H,
                    n = nrow(galaxy),
                    y = galaxy$k_speed,
                    grid_max = grid_max,
                    grid_min = grid_min,
                    grid_length = grid_length),
                iter = 2000,
                chains = 12*3,
                seed = 900817669)
## Warning: There were 272 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 5167 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Check lp__ traces to detect better fits
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::gather_draws(lp__) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Order fits by mean lp__ to find bad ones
normal_fixed_mixture_stan_tau6_good <-
    normal_fixed_mixture_stan_tau6 %>%
    tidybayes::gather_draws(lp__) %>%
    group_by(.chain) %>%
    summarize(mean = mean(.value)) %>%
    arrange(mean) %>%
    ## Define
    mutate(good_chain = between(mean, max(mean) - 1, max(mean)))
normal_fixed_mixture_stan_tau6_good %>%
    print(n = Inf)
## # A tibble: 36 x 3
##    .chain  mean good_chain
##     <int> <dbl> <lgl>     
##  1      7 -217. FALSE     
##  2      1 -217. FALSE     
##  3     20 -217. FALSE     
##  4     17 -217. FALSE     
##  5     26 -217. FALSE     
##  6     16 -217. FALSE     
##  7     30 -217. FALSE     
##  8     31 -217. FALSE     
##  9     14 -217. FALSE     
## 10     10 -217. FALSE     
## 11     24 -217. FALSE     
## 12      2 -217. FALSE     
## 13     34 -216. FALSE     
## 14      5 -216. FALSE     
## 15     21 -216. FALSE     
## 16      6 -216. FALSE     
## 17      3 -216. FALSE     
## 18      9 -216. FALSE     
## 19     32 -216. FALSE     
## 20     13 -216. FALSE     
## 21     12 -216. FALSE     
## 22     28 -216. FALSE     
## 23      4 -216. FALSE     
## 24     33 -216. FALSE     
## 25     19 -216. FALSE     
## 26     18 -216. FALSE     
## 27     22 -216. FALSE     
## 28     36 -216. FALSE     
## 29     27 -216. TRUE      
## 30     25 -216. TRUE      
## 31     11 -216. TRUE      
## 32     35 -215. TRUE      
## 33     23 -215. TRUE      
## 34     15 -215. TRUE      
## 35     29 -215. TRUE      
## 36      8 -215. TRUE
## For vague prior fit
bad_chains <- normal_fixed_mixture_stan_tau6_good %>%
    filter(!good_chain) %>%
    magrittr::extract2(".chain")

## Check lp__ traces for good chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::gather_draws(lp__) %>%
    filter(!(.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Check lp__ traces for bad chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::gather_draws(lp__) %>%
    filter(.chain %in% bad_chains) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Rhat for good chains (slightly different definition?)
lapply(normal_fixed_mixture_stan_tau6@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,-1 * bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 1.002634
## Rhat among bad chains
lapply(normal_fixed_mixture_stan_tau6@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 1.006507
## Estimated densities among good chains
plot_draws(normal_fixed_mixture_stan_tau6,
           chains = setdiff(seq_len(length(normal_fixed_mixture_stan_tau6@sim$samples)),
                            bad_chains))

plot_draws(normal_fixed_mixture_stan_tau6,
           chains = setdiff(seq_len(length(normal_fixed_mixture_stan_tau6@sim$samples)),
                            bad_chains)[1])

## Scatter plot of parameter values
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::tidy_draws() %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = `mu[1]`, y = `mu[2]`)) +
    geom_point() +
    facet_grid(~ bad_chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among good chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among bad chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among good chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among bad chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among good chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among bad chains
normal_fixed_mixture_stan_tau6 %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

pairs_plot_all(normal_fixed_mixture_stan_tau6)

pairs(normal_fixed_mixture_stan_tau6, pars = c("mu","sigma","Pi"))

plot_draws(normal_fixed_mixture_stan_tau6)

## shinystan::launch_shinystan(normal_fixed_mixture_stan_tau6)

Finite Mixture of Normals (Ordered mu)

3 latent cluster model

H <- 3
normal_fixed_mixture_stan_tau3_order <-
    rstan::stan(model_code = readr::read_file("./bayesianideas_density_normal_fixed_mixture_tau.stan"),
                data = list(
                    ## tau ~ Gamma(alpha, beta)
                    alpha = rep(10^(-3), H), beta = rep(10^(-3), H),
                    ## mu ~ N(m, s_squared)
                    m = rep(20,H), s_squared = rep(1000,H),
                    ## Pi ~ Dirichlet(dirichlet_alpha / H)
                    dirichlet_alpha = rep(1,H),
                    ## Number of clusters
                    H = H,
                    n = nrow(galaxy),
                    y = galaxy$k_speed,
                    grid_max = grid_max,
                    grid_min = grid_min,
                    grid_length = grid_length),
                iter = 2000,
                chains = 12*3,
                seed = 339492111)
## Warning: There were 191 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 70 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Check lp__ traces to detect better fits
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::gather_draws(lp__) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Order fits by mean lp__ to find bad ones
normal_fixed_mixture_stan_tau3_order_good <-
    normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::gather_draws(lp__) %>%
    group_by(.chain) %>%
    summarize(mean = mean(.value)) %>%
    arrange(mean) %>%
    ## Define
    mutate(good_chain = between(mean, max(mean) - 1, max(mean)))
normal_fixed_mixture_stan_tau3_order_good %>%
    print(n = Inf)
## # A tibble: 36 x 3
##    .chain  mean good_chain
##     <int> <dbl> <lgl>     
##  1      1 -225. FALSE     
##  2     25 -223. FALSE     
##  3     27 -222. FALSE     
##  4     13 -222. FALSE     
##  5     20 -222. FALSE     
##  6      3 -221. FALSE     
##  7     24 -221. FALSE     
##  8      7 -220. FALSE     
##  9     28 -218. FALSE     
## 10     15 -216. FALSE     
## 11     23 -215. FALSE     
## 12     29 -215. FALSE     
## 13     19 -215. FALSE     
## 14     32 -215. FALSE     
## 15     21 -215. FALSE     
## 16     36 -213. FALSE     
## 17     33 -213. FALSE     
## 18      4 -213. FALSE     
## 19     16 -213. FALSE     
## 20     26 -213. FALSE     
## 21      5 -213. FALSE     
## 22     35 -213. FALSE     
## 23      6 -212. FALSE     
## 24     14 -206. TRUE      
## 25     31 -206. TRUE      
## 26     10 -206. TRUE      
## 27     17 -206. TRUE      
## 28      9 -206. TRUE      
## 29     34 -206. TRUE      
## 30     12 -206. TRUE      
## 31      2 -206. TRUE      
## 32     22 -206. TRUE      
## 33     30 -206. TRUE      
## 34      8 -206. TRUE      
## 35     18 -206. TRUE      
## 36     11 -206. TRUE
## For vague prior fit
bad_chains <- normal_fixed_mixture_stan_tau3_order_good %>%
    filter(!good_chain) %>%
    magrittr::extract2(".chain")

## Check lp__ traces for good chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::gather_draws(lp__) %>%
    filter(!(.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Check lp__ traces for bad chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::gather_draws(lp__) %>%
    filter(.chain %in% bad_chains) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Rhat for good chains (slightly different definition?)
lapply(normal_fixed_mixture_stan_tau3_order@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,-1 * bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 1.002835
## Rhat among bad chains
lapply(normal_fixed_mixture_stan_tau3_order@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 1.643682
## Estimated densities among good chains
plot_draws(normal_fixed_mixture_stan_tau3_order,
           chains = setdiff(seq_len(length(normal_fixed_mixture_stan_tau3_order@sim$samples)),
                            bad_chains))

## Scatter plot of parameter values
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::tidy_draws() %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = `mu[1]`, y = `mu[2]`)) +
    geom_point() +
    facet_grid(~ bad_chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among good chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among bad chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among good chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among bad chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among good chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among bad chains
normal_fixed_mixture_stan_tau3_order %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

pairs_plot_all(normal_fixed_mixture_stan_tau3_order)

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

pairs(normal_fixed_mixture_stan_tau3_order, pars = c("mu","sigma","Pi"))
## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

plot_draws(normal_fixed_mixture_stan_tau3_order)

## shinystan::launch_shinystan(normal_fixed_mixture_stan_tau3_order)

6 latent cluster model

H <- 6
normal_fixed_mixture_stan_tau6_order <-
    rstan::stan(model_code = readr::read_file("./bayesianideas_density_normal_fixed_mixture_tau.stan"),
                data = list(
                    ## tau ~ Gamma(alpha, beta)
                    alpha = rep(10^(-3), H), beta = rep(10^(-3), H),
                    ## mu ~ N(m, s_squared)
                    m = rep(20,H), s_squared = rep(1000,H),
                    ## Pi ~ Dirichlet(dirichlet_alpha / H)
                    dirichlet_alpha = rep(1,H),
                    ## Number of clusters
                    H = H,
                    n = nrow(galaxy),
                    y = galaxy$k_speed,
                    grid_max = grid_max,
                    grid_min = grid_min,
                    grid_length = grid_length),
                iter = 2000,
                chains = 12*3,
                seed = 410921792)
## Warning: There were 423 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 18410 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Check lp__ traces to detect better fits
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::gather_draws(lp__) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Order fits by mean lp__ to find bad ones
normal_fixed_mixture_stan_tau6_order_good <-
    normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::gather_draws(lp__) %>%
    group_by(.chain) %>%
    summarize(mean = mean(.value)) %>%
    arrange(mean) %>%
    ## Define
    mutate(good_chain = between(mean, max(mean) - 1, max(mean)))
normal_fixed_mixture_stan_tau6_order_good %>%
    print(n = Inf)
## # A tibble: 36 x 3
##    .chain  mean good_chain
##     <int> <dbl> <lgl>     
##  1     19 -214. FALSE     
##  2     11 -208. FALSE     
##  3     26 -208. FALSE     
##  4      3 -208. FALSE     
##  5      6 -207. FALSE     
##  6     12 -207. FALSE     
##  7     27 -207. FALSE     
##  8      5 -207. FALSE     
##  9     35 -207. FALSE     
## 10     18 -207. FALSE     
## 11     16 -207. FALSE     
## 12     29 -207. FALSE     
## 13     22 -207. FALSE     
## 14      7 -207. FALSE     
## 15     15 -207. FALSE     
## 16     34 -207. FALSE     
## 17     25 -207. FALSE     
## 18     10 -206. FALSE     
## 19     30 -206. FALSE     
## 20     36 -206. FALSE     
## 21      1 -206. FALSE     
## 22     13 -206. FALSE     
## 23     14 -206. FALSE     
## 24      8 -206. FALSE     
## 25     23 -206. FALSE     
## 26      9 -206. FALSE     
## 27     20 -206. FALSE     
## 28     21 -206. FALSE     
## 29     24 -206. FALSE     
## 30      2 -206. FALSE     
## 31     33 -206. FALSE     
## 32     32 -206. FALSE     
## 33     28 -206. FALSE     
## 34     17 -205. TRUE      
## 35      4 -204. TRUE      
## 36     31 -204. TRUE
## For vague prior fit
bad_chains <- normal_fixed_mixture_stan_tau6_order_good %>%
    filter(!good_chain) %>%
    magrittr::extract2(".chain")

## Check lp__ traces for good chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::gather_draws(lp__) %>%
    filter(!(.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Check lp__ traces for bad chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::gather_draws(lp__) %>%
    filter(.chain %in% bad_chains) %>%
    ggplot(mapping = aes(x = .iteration, y = .value, group = .chain)) +
    geom_line() +
    facet_wrap( ~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## Rhat for good chains (slightly different definition?)
lapply(normal_fixed_mixture_stan_tau6_order@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,-1 * bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 1.024318
## Rhat among bad chains
lapply(normal_fixed_mixture_stan_tau6_order@sim$samples,
       "[[", "lp__") %>%
    do.call(cbind, .) %>%
    "["(,bad_chains) %>%
    asbio::R.hat(M = ., burn.in = 0.5)
## [1] 1.082194
## Estimated densities among good chains
plot_draws(normal_fixed_mixture_stan_tau6_order,
           chains = setdiff(seq_len(length(normal_fixed_mixture_stan_tau6_order@sim$samples)),
                            bad_chains))

## Scatter plot of parameter values
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::tidy_draws() %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    ggplot(mapping = aes(x = `mu[1]`, y = `mu[2]`)) +
    geom_point() +
    facet_grid(~ bad_chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among good chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## mean traces among bad chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::spread_draws(mu[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = mu, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among good chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## precision traces among bad chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::spread_draws(tau[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = tau, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among good chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(!bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

## proportion traces among bad chains
normal_fixed_mixture_stan_tau6_order %>%
    tidybayes::spread_draws(Pi[h]) %>%
    mutate(bad_chain = (.chain %in% bad_chains)) %>%
    filter(bad_chain) %>%
    ungroup() %>%
    ggplot(mapping = aes(x = .iteration, y = Pi, group = .chain, color = factor(h))) +
    geom_line() +
    facet_wrap(~ .chain) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

pairs_plot_all(normal_fixed_mixture_stan_tau6_order)

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

pairs(normal_fixed_mixture_stan_tau6_order, pars = c("mu","sigma","Pi"))
## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, : Binning grid too coarse for current
## (small) bandwidth: consider increasing 'gridsize'

plot_draws(normal_fixed_mixture_stan_tau6_order)

## shinystan::launch_shinystan(normal_fixed_mixture_stan_tau6_order)
print(sessionInfo())
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.4
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bayesplot_1.6.0    tidybayes_1.0.3    DPpackage_1.1-7.4  rstan_2.18.2       StanHeaders_2.18.1
##  [6] forcats_0.3.0      stringr_1.4.0      dplyr_0.8.0        purrr_0.3.0        readr_1.3.1       
## [11] tidyr_0.8.2        tibble_2.0.1       ggplot2_3.1.0      tidyverse_1.2.1    doRNG_1.7.1       
## [16] rngtools_1.3.1     pkgmaker_0.27      registry_0.5       doParallel_1.0.14  iterators_1.0.10  
## [21] foreach_1.4.4      knitr_1.21        
## 
## loaded via a namespace (and not attached):
##  [1] pixmap_0.4-11             nlme_3.1-137              matrixStats_0.54.0        lubridate_1.7.4          
##  [5] httr_1.4.0                tools_3.5.2               backports_1.1.3           utf8_1.1.4               
##  [9] R6_2.4.0                  KernSmooth_2.23-15        lazyeval_0.2.1            colorspace_1.4-0         
## [13] withr_2.1.2               tidyselect_0.2.5          gridExtra_2.3             prettyunits_1.0.2        
## [17] processx_3.2.1            compiler_3.5.2            cli_1.0.1                 rvest_0.3.2              
## [21] asbio_1.5-5               arrayhelpers_1.0-20160527 xml2_1.2.0                labeling_0.3             
## [25] scales_1.0.0              mvtnorm_1.0-8             ggridges_0.5.1            callr_3.1.1              
## [29] multcompView_0.1-7        digest_0.6.18             rmarkdown_1.11            pkgconfig_2.0.2          
## [33] htmltools_0.3.6           bibtex_0.4.2              plotrix_3.7-4             rlang_0.3.1              
## [37] readxl_1.2.0              rstudioapi_0.9.0          generics_0.0.2            svUnit_0.7-12            
## [41] jsonlite_1.6              inline_0.3.15             magrittr_1.5              loo_2.0.0                
## [45] Matrix_1.2-15             Rcpp_1.0.0                munsell_0.5.0             fansi_0.4.0              
## [49] scatterplot3d_0.3-41      stringi_1.3.1             yaml_2.2.0                MASS_7.3-51.1            
## [53] pkgbuild_1.0.2            plyr_1.8.4                ggstance_0.3.1            grid_3.5.2               
## [57] crayon_1.3.4              lattice_0.20-38           haven_2.0.0               splines_3.5.2            
## [61] hms_0.4.2                 ps_1.3.0                  pillar_1.3.1              tcltk_3.5.2              
## [65] reshape2_1.4.3            codetools_0.2-16          stats4_3.5.2              glue_1.3.0               
## [69] evaluate_0.13             modelr_0.1.3              deSolve_1.21              cellranger_1.1.0         
## [73] gtable_0.2.0              assertthat_0.2.0          xfun_0.4                  xtable_1.8-3             
## [77] broom_0.5.1               coda_0.19-2               survival_2.43-3
## Record execution time and multicore use
end_time <- Sys.time()
diff_time <- difftime(end_time, start_time, units = "auto")
cat("Started  ", as.character(start_time), "\n",
    "Finished ", as.character(end_time), "\n",
    "Time difference of ", diff_time, " ", attr(diff_time, "units"), "\n",
    "Used ", foreach::getDoParWorkers(), " cores\n",
    "Used ", foreach::getDoParName(), " as backend\n",
    sep = "")
## Started  2019-05-11 09:38:18
## Finished 2019-05-11 10:12:10
## Time difference of 33.86934 mins
## Used 12 cores
## Used doParallelMC as backend